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Abstract 

Using variational density matrix optimization with two- and three-index conditions we study 
the one-dimensional Hubbard model with periodic boundary conditions at various filling factors. 
Special attention is directed to the full exploitation of the available symmetries, more specifi- 
cally the combination of translational invariance and space-inversion parity, which allows for the 
study of large lattice sizes. We compare the computational scaling of three different semidefinite 
programming algorithms with increasing lattice size, and find the boundary point method to be 
the most suited for this type of problem. Several physical properties, such as the two-particle 
correlation functions, are extracted to check the physical content of the variationally determined 
density matrix. It is found that the three-index conditions are needed to correctly describe the full 
phase diagram of the Hubbard model. We also show that even in the case of half filling, where 
the ground-state energy is close to the exact value, other properties such as the spin-correlation 
function can be flawed. 
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INTRODUCTION 



The reduced density matrix makes its first appearance in the work of Dirac, in which 
the single-particle density matrix (IDM) is used in the description of Hartree-Fock theory 
[1]. Husimi |2] was the first to note that, for a system of identical particles interacting only 
in a pairwise manner, the energy can be expressed exactly as a function of the 2DM. This 
becomes very clear in second-quantized notation (see e.g. [3[|3]), where a system of identical 
particles interacting pairwise is described by a general Hamiltonian: 

= ^ tapaiap + yai3;-f5oialasa^ . (1) 

a/3 afS-fS 

The expectation value for the energy of any ensemble of A^-particle wave functions {"^f) 
with positive weights Wi, can then be expressed as a function of the 2DM alone: 

Y: w^i^mm = Tr ri7(^) = \ r./3;7.^aW ' (2) 

i af3-fS 

in which we have introduced the 2DM: 

ra/3;75 = y^^Wi{'^^\aialasa-y\'^^) , with ^Wj = l, (3) 

i i 

and the reduced two-particle Hamiltonian, 

(2) 

-^a/3;7(5 ~ ]\[ — I ^^^^^1^^ ~ ^aStp-i — SfSjtaS + Spsta-f) + Va/Si-yS ■ (4) 

The idea to use the 2DM as a variable in a variational scheme was first published in literature 
by Lowdin in his groundbreaking article [3], but even earlier, in 1951, John Coleman tried 
a practical variational calculation on Lithium. To his surprise, the energy he obtained was 
far too low, after which he realized the variation was performed over too large a class of 
2DM's [6]. Independently and unaware of the work by Lowdin and Coleman, Joseph Mayer 
[7] used the 2DM in a study of the electron gas. In a reply to Mayer's paper, Tredgold [8] 
pointed out the unphysical nature of the results, and suggested that additional conditions 
on the density matrix are needed to improve on them. 

These results led Coleman, in his seminal review paper [9J, to formulate the A^-represen- 
tability problem. This is the problem of finding the necessary and sufficient conditions which 
a reduced density matrix has to fulfil to be derivable from a statistical ensemble of physical 
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wave functions, i. e. expressible as in Eq. In this paper he also derived the necessary and 
sufficient conditions for ensemble iV-represent ability of the IDM, and some bounds on the 
eigenvalues of the 2DM. A big step forward was the derivation of the Q and Q matrix non- 
negativity conditions by Garrod and Percus [lOj. These were practical constraints, which 
allowed for a computational treatment of the problem. The ffist numerical calculation using 
these conditions on the Beryllium atom [HI [12] was very encouraging, as the results were 
highly accurate. It turned out, however, that Beryllium, due to its simple electronic stucture, 
is a special case where these conditions perform very well. A subsequent study showed that 
these conditions do not work well at all for nuclei [T3|[Ti]. This disappointing result, together 
with the computational complexity of the problem, caused activity in the field to diminish 
for the next 25 years. The change came with the development of a new numerical technique, 
called semidefinite programming, which turned out to be very suited for the determination 
of the 2DM under matrix non-negativity constraints. Maho Nakata et al. [15] were the ffist 
to use a standard semidefinite programming package to calculate the ground-state energies 
of some atoms and molecules, and obtained quite accurate results. He was quickly followed 
by the extensive work of Mazziotti [16j. These results reinvigorated interest in the method, 
and sparked of a lot of developments. New A^-represent ability conditions were introduced, 
e.g. the three-index T conditions, set forth by Zhao et al. [17], which led to mHartree 
accuracy [T81 - I23] for molecules near equilibrium geometries. 

In recent years interest in the method has been growing, as the variational determination 
of the 2DM results in a lower bound, which is highly complementary to the upper bound 
obtained in variational approaches based on a wave-function ansatz. In addition, the method 
is essentially non-perturbative in nature, and has a completely different structure quite 
unrelated to other many-body techniques. A lot of activity has been devoted to the search for 
new A^-representability conditions, which improve the result in a computationally cheap way 
[211 ES] ■ There have also been efforts to improve the semidefinite programming algorithms 
by adapting them to the specific problem of density matrix optimization [261428] . allowing 
the study larger systems. 

The one-dimensional Hubbard model p9] is the simplest model possessing non-trivial 
correlations present in a solid. The Hamiltonian reads: 

ia i 
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It pertains to a one-dimensional lattice, with sites labeled by i = 1, . . . , L. Periodic boundary 
conditions (PBC) are assumed. 

The complexity of this seemingly simple schematic Hamiltonian lies in the competition 
between the first term, called the hopping term, which delocalizes the electrons, and the 
second on-site repulsion term which is diagonal in the site basis. In Figure [T] a graphic 
representation of the two terms in the Hubbard Hamiltonian is shown. 

In this article the one-dimensional Hubbard is studied using the standard XQQ two-index 
and 7i,2 three-index non-negativity conditions. In the next Section we discuss how the huge 
amount of symmetry which is present in the model can be exploited to get a significant 
speedup of the programs. More specifically, it is shown how translational invariance and 
space-inversion parity can be combined. In the subsequent Section this symmetry is used to 
compare the computational scaling of different semidefinite programming algorithms with 
increasing lattice sites. In the final Section we show the v2DM results using two- and 
three-index constraints for the one-dimensional Hubbard model on a 20- and 50-site lattice 
with different filling factors. We have not only computed the ground-state energy, but 
also compared the spin and charge two-particle correlations functions with Quantum Monte 
Carlo [30] and Bethe ansatz |31] results, to check the validity of the variationally obtained 
2DM. 

EXPLOITING SYMMETRY IN THE HUBBARD MODEL 

As there is no preferred direction in spin space in Eq. (|5]), spin symmetry can be exploited 
and all the results from [32] are taken over. However, far more symmetries are present in 
the model, which allow for a huge reduction of the matrix dimensions involved. 

Translational invariance 

When periodic boundary conditions are assumed (as in Figure [T]) the Hamiltonian is 
invariant under translations along the lattice. This is an abelian symmetry for which it is 
easy to block diagonalize the 2DM and its matrix maps, as the correct basis transformation 
in single-particle space automatically block diagonalizes all matrices on higher-order particle 
space. Translational invariance is exploited by Fourier transforming the site basis to quasi- 
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FIG. 1: Illustration of the two terms present in the Hubbard Hamiltonian, electrons can jump to 
nearest-neighbour sites with amplitude tij. When two electrons are on the same site, there is an 
energy penalty of U. 

momentum space: 

1^^) = VzE^^''!-^'^) , (6) 

j 

where L is the lattice size, and k takes on the values ^ for n = 0, . . . , L — 1. The kinetic 
or hopping part of the Hamiltonian becomes diagonal in this basis: 

-^hop = -2t COS A; al^ttka , (7) 

ka 

from which it follows that in the non-interacting {U = 0) ground state, the electrons occupy 
the states with momenta lower than the Fermi level. 

The eigenstates of the Hubbard Hamiltonian have a good total quasi- momentum /C. The 
2DM for these states, expressed in the quasi-momentum single-particle basis: 

^kakh;kckd = '^'^iTc]2'^^'^SM,i\B\akt ^kcka\'^SM,i) ' i^) 

5 



where 



V cr;,(T(, 



is automatically block diagonal, because the only non-zero matrix elements in Eq. (|8j) are 
those which conserve momentum: {ka + /c;,)%27r = {kc + kd)%27r^^. This means we have L 
blocks T^^ for every 5*, with two-particle states that satisfy K = {ka + kh)%27T, and a block 
dimension that scales linearly with lattice size L. 

The spin-symmetric matrix constraints simplify considerably by including translational 
invariance, because the IDM is automatically diagonal in the quasi-momentum basis: 



Pk 



(10) 



M 



The translationally invariant IDM can be derived from the 2DM as: 

P. = ^ E 5^^' + 4.0rf^.,, . (11) 

5 k' 

As an example the translationally invariant form of the Q condition is shown. There is a 
slight complication because the correct annihilation or hole operator is given by: 



O-ka 



:-l)i+'^a 



k—a ; 



with k 



327r . 



(12) 



Using this operator the translationally invariant Q map becomes: 



SK 



^^12 2-^\^SM,i\^ kakt ^kcka\^SM,i/ 



(13) 



where K = {ka + kb)%2n = {kc + fc(i)%27r and with the particle-hole operator defined by: 



At 



SK 

kn kh 



SK 



(14) 



The Q map can be expressed as a function of the 2DM: 



1 1 



S 



Si^)kM;kcka - KkAak.Pka - Y (1 + ^ka~ki){^ + ^k„k,) \ \ [ ^kah;kch ' (^^) 



from which one can see that the blocks in the Q matrix with K = {ka + kh)%27T = {kc + 
fcd)%27r, correspond to the blocks with K' = {ka + kd)%2TT = {kc + kb)%2T^ in the 2DM, as 
expected for a particle-hole transformed quantity. 
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Parity 



The Hubbard model with periodic boundary conditions (PBC) has more symmetries than 
spin and translational invariance, one of them being parity. Parity follows from the symmetry 
under the inversion of space, i.e. r — t- — r. One can readily appreciate from Figure [T] that 
this symmetry is present here, i.e. the Hamiltonian is invariant for the inversion of the site 
index i — t- —i%L. From the Fourier transform in Eq. ^ it can be seen that the effect of this 
operator on a momentum state is to transform k into k = —k%2TT. However, this operation 
does not commute with translation, which means that the eigenstates of the Hamiltonian 
cannot have good momentum and parity at the same time. From now on we assume that 
the number of lattice sites L is even. 

As a consequence, if the ground state has momentum |/C) 7^ or tt then it is doubly 
degenerate, forming a doublet with |/C). In what follows we will use this degeneracy to 
exploit both translational invariance and parity to reduce the dimensions of the matrices 
involved in the program. The following considerations are valid for every /C. We start with 
the simplest case, the IDM. 

Single-particle space is built out of L different momentum states having up or down spin, 
\ka), for < k < 271 and a = ±|. If we transform to a basis with good parity, momentum 
is no longer a good quantum number: 

|A;^(t) = PNk {\ka) + n\k(T)) , with < < tt . (16) 

Two states, k = and k = it [3^] are mapped on themselves, and only positive parity 
states can be formed with these momenta. They have norms ^N^ = \- For the other states, 
< A; < TT, both positive and negative parity combinations can be constructed, with norm 
PNk = To take advantage of both symmetries at the same time, we define the IDM 
using an ensemble of the |/C) and |^) states: 

Pu^l'^' =E^4l^5:E(^S5?l4.a^-.I^S?) (17) 



ml 



2 ^ 



E [^^^' [pf + ^^'/^f ) + + ^/^f ) ] ' (18) 



in which ffj^ is a regular translationally invariant IDM as defined in Eq. (10). Because of 
parity symmetry we have. 



Pk=ph (19) 



from which it follows that Eq. (18) is diagonal in vr. If we now define: 
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Pk = ^ [p^ + pI) , (20) 
the translationally invariant IDM with good parity reduces to: 

Pk- = \{pk + Pi) , (21) 

which can be seen to be independent of parity. The IDM is still diagonal in k, but less 
elements have to be stored, since for < A; < tt there is a degeneracy in parity. 

In contrast to parity in atomic systems (see e.g. L32j), the parity of a two-particle state 
for translationally invariant systems is not the product of the parities of the single-particle 
states building up the two-particle state. Instead, the parity is inherently a two-particle 
property, and the single-particle states building up the two-particle states have no good 
parity: 

\ab; Sk"") = ^N^^ {\ab; SK) + 7T\ab; SK)) , (22) 

with 

< ^ < TT and < a, 6 < 27r . (23) 
In general the 2DM is defined using a pK, ensemble: 

'-ab^cd- Z^^i 2^ l^^^SMA^-^ab ^cd \^ SM,i I ^ (24) 

where 

B'T = 'N^b [b^T + -B^f) , (25) 

and with B"^ defined as in Eq. ( 9 ) . Because of this p/C-ensemble definition and the fact 
that parity symmetry implies that, 

IC-pSK _ -IC-pSK _ /9R\ 
J- ab;cd — ab;cd ' K"^^) 

one sees that the 2DM becomes diagonal in two-particle parity. As was the case for the 
IDM, the K = and vr are mapped on themselves, but because the single-particle momenta 
a, b change, both positive and negative parity combinations can now be formed. Let us take 
a look at the different possibilities: 
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K = 0; ioT K = the single-particle indices in Eq. (22) have to satisfy: 

(a + b)%2n = or a = b . (27) 
This means that K = states can be written as: 



\aa 



■S0'')=^Nj^^{\aa;S0) + 7c\aa;S0)) , (28) 



in which the second term is equal to the first but with exchanged single-particle indices. 
From previous discussions we know that the symmetry of the two-particle state under the 
exchange of the single-particle indices depends on the two-particle spin, i. e. 



aa;SO) = (-l)'^|aa; ^0) . (29) 



One can see from Eq. (28) that for 5* = only the positive parity states, and for S = 1 
only negative parity states remain. The norm is given ^N^^, ~ K = case, the 

definition of the parity-symmetric 2DM as a function of the regular translationally invariant 
2DM then reduces to: 

rf°;cd = '^7r(-i)srf°.c^ . (30) 



< K < 77 .■ for < K < 71 the first and second term in Eq. ( 22 ) consist of different 
single-particle indices a y^b, implying that both positive and negative parity combinations 
can be constructed, with norm '"A'"^ = As shown for the IDM, the pJC ensemble makes 
the 2DM diagonal in, and independent of, parity. Hence every block is twofold degenerate. 
Since K ^ K there are only two terms remaining in the definition of the parity-symmetric 
2DM: 

ab;cd ~ 2 V "■^'^'^ ) ' ^ 

K = 77.- Finally, for this block K again equals K. In this case there is always one state 
that is mapped on itself, and for which only a positive parity combination can be formed, 
i.e. |07r;S'7r+), with norm '"A*",^^ = |. For all the other states in this block both positive 
and negative parity combinations can be formed, with norm ^N^i, = Because of the p/C 
ensemble, the 2DM falls apart in a positive and negative parity block, and since K = K, 
four terms remain in the definition of the 2DM: 

J- ab;cd — ^^ab cd y'- ab;cd ^ '- ab;cd ^ ^ [J- ab;cd ^ ^ ab;cd\ J ' K"^"^) 
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We observe from Eq. (32) that the original K = n block reduces to a positive and negative 
parity block, for both the S = and S = 1 part. Also note that there is no degeneracy 
between the positive and negative parity block! 

Similar considerations hold for the matrix constraints, as an example the explicit case of 
the Q condition is given. 

The parity-symmetric form of a particle-hole state is defined as: 

\ab;Sk^) = '^N^, [ya®~a,f''+[ai®a-i!\''^ |0) , (33) 



in which the hole operator ata is defined as in Eq. ( 12 ). Using this parity-symmetric particle- 
hole operator the Q map is defined in a p/C ensemble, which once again renders the matrix 
diagonal in particle-hole parity. The particle-hole states can be divided into two classes, 
on the one hand = or tt, and on the other hand those states which are mapped on a 
different momentum. For simplicity we first consider this last class. 

< K < 77.- in this case one can construct both positive and negative parity combina- 
tions, with norm ^A*"^ = The resulting Q matrix contains only two terms, because of 
momentum conservion, and is independent of particle-hole parity, so every block is twofold 
degenerate: 
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^(r)^^;, = - [^(r)^^,, + G{T)%,\ . (34) 

This implies the following expression of the Q map as a function of the parity-symmetric 
2DM: 

y in ab;cd=^acdbdpa 4 FpfK' Tj^K' 1^^^ \ \^ , ^, ( 1^^ ad;cb ■ (35) 

ad cb S' 2 2 J t' 

K = and K = tt.- both the K = and K = tt blocks are mapped on themselves. For 
K = the action of the parity operator is again to exchange the single-particle momenta, 
but in contrast with the two-particle case, there is no symmetry between the particle and 
the hole index. As a consequence positive and negative parity combinations for both K = 
and K = TT can be constructed, with norms = There are a few exceptions however: 
for A' = 0, the states with a = 6 = and a = b = tt, and for K = n the states with a = 0, 
b = TT and a = tt, 6 = 0, are mapped on themselves and only occur in the positive parity 
block, with norm = ^. The general form of the parity-symmetric Q map when K = K, 
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FIG. 2: The number of basic matrix operations needed to converge for different lattice sizes of 
the one-dimensional Hubbard model, using XQQ conditions and U = 1. 



as a function of the regular translationally invariant Q is: 



)SK 
ab;cd 



(36) 



In this case the expression of ^ as a function of the 2DM is a bit more complicated: 



1 1 



s 



11c/ 
2 2'-' 



V(l + ^ad)(l + '^e5) yS'l 



cb 



ad cb 



+7r 



^{l + 6ad){l + S,t) 



4 'N^' 'Kb 



-r ucb) S^-pS' 

r ]\fK" ^ ad 



dd\cb 



(37) 



COMPUTATIONAL PERFORMANCE 



In v2DM we want to optimize the energy by varying a matrix, the 2DM, under the 
constraints that it has the right particle number, and that some linear matrix maps of the 
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2DM are positive semi definite, i.e. 



ininTr \rH^^^] 

( N(N-1) 




(38) 



u.c.t. 



2 



(39) 



A- (r) ^ . 



Here the Cj are a collection of matrix non-negativity constraints to be applied. It is well 
known that this problem can be reformulated as a semidefinite program. There is a vast 
literature on this subject and many different algorithms exist to solve this type of problem. 
Because of the large amount of symmetry present in the Hubbard model, there is a huge block 
diagonalization of the 2DM and the matrix constraints. This leads to a significant reduction 
of the computational cost of a basic matrix computation. This feature has allowed us to go 
up to large lattice sizes and compare the computational scaling of different algorithms. We 
have implemented three different semidefinite programming algorithms. Two are so-called 
interior point methods (a dual-only potential reduction algorithm, [32] and a primal-dual 
interior point algorithm [28j), where the 2DM is optimized from within the A^-representable 
region. In the third one (a boundary point method p7j) the 2DM is not required to be 
iV-representable during the optimization. In this Section we compare the computational 
scaling of these methods for the one-dimensional Hubbard model with U = 1 and using 
XQQ conditions. For details about the implementation of the different algorithms, see the 
cited references. 

All the methods have the same basic computational scaling behaviour, being O(M^) 
(with M the dimension of single-particle Hilbert space) required for multiplying, inverting 
or diagonalizing a matrix. In Figure |2] the number of these operations needed to converge to 
the optimum is plotted as a function of lattice size. The interior point methods both have to 
solve a linear system of size M^, so it is not surprising that the scaling, on top of the basic 
matrix computations, of these methods is M^. More surprising is that there seems to be no, 
or a very limited, scaling for the boundary point method. The number of iterations required 
remains around 3000, irrespective of the size of the system. It must be stressed that this is 
a result limited to the one-dimensional hubbard model, and cannot be extrapolated to other 
systems, like molecules, where the convergence properties of the boundary point method can 
be completely different. One reason for the succes of the boundary point method applied to 
the Hubbard model is the amount of symmetry present in the system. The boundary point 
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method is designed for problems with a huge amount of dual variables or primal constraints. 
For most physical systems the dimensions of the matrices involved are already unfeasibly 
large before the boundary point method would becomes advantageous. The Hubbard model, 
however, contains many symmetries, implying that the matrix dimensions are considerably 
reduced, and the number of dual variables can get very large before the matrix computations 
involved become unfeasible. In this case, the domain where the boundary point method is 
advantageous is actually reached. 

RESULTS 

In this Section we present and discuss the results of v2DM calculations, taking advan- 
tage of all the symmetries, on a 50-site lattice with the XQQ conditions, and on a 20-site 
lattice with the XQQT1T2 (XQQT) conditions. The Hubbard model has been studied be- 
fore using the v2DM method, see e.g. [191 I33ti35] . but only the energy was considered, 
and this for relatively small lattice sizes (up to L = 14). In this paper we study different 
filling factors, and extract various properties like the ground-state energy and two-particle 
correlation functions in order to assess the quality of the variationally obtained 2DM. The 
v2DM results discussed in this Section were all obtained using the primal-dual predictor 
corrector semidefinite programming algorithm [2B] . Although the one-dimensional Hubbard 
model can be solved exactly using the Bethe ansatz [36ti39] , it is hard to extract information 
about the solution for finite systems. For the calculations on a 20-site lattice, we compare 
the data with the quasi-exact results obtained through a variational Matrix Product State 
(MPS) algorithm |^0] - H2] . written by co-worker Sebastian Wouters [13]. For the 50-site lat- 
tice, however, this is no longer computationally feasible. At half filling a simplification in the 
Bethe-ansatz equations occurs, which allows to calculate the ground-state energy of finite 
systems by solving a set of non-linear equations (Lieb-Wu)[44]. At other fillings no data is 
available for comparison. 

Ground-state energy 

In Fig. [3] the ground-state energy per particle of the one-dimensional Hubbard model 
is plotted as a function of the on-site repulsion U (the hopping parameter t will always 
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FIG. 3: Ground-state energy per particle as a function of on-site repulsion U of the Hubbard model 
for a 20-site (top) and 50-site (bottom) lattice at half-, ^ and ^ filling. For the 20-site lattice a 
comparison is made between v2DM using the TQQ and XQQT conditions, and a quasi-exact result 
using MPS. For the 50-site lattice only XQQ conditions are feasible, and these have been compared 
to the exact (Bethe-ansatz) result for half filling. 
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be taken equal to unity). In the top figure the v2DM results for the 20-site lattice are 
shown for three different fillings, 12 particles (^), 16 particles (^) and half filling. These 
were calculated using both the XQQ and the XQQT conditions, and are compared to the 
quasi-exact variational MPS results. In the bottom figure the v2DM results for the 50- 
site lattice are shown for the same fillings {i.e. 30 particles {j^), 40 particles and half 
filling). For the 50-site lattice it was only possible to perform the calculations using the XQQ 
conditions, and compare to the exact solution obtained by solving the Lieb-Wu equations 
for the half-filled lattice [H]. 

One interesting thing to notice is that the XQQ energy per particle for the 20-site lattice 
and the 50-site lattice, at the same filling, are very similar. This is due to the periodic 
boundary conditions which make the results converge quite rapidly for increasing lattice 
size L, implying that one can already extract relevant results for the thermodynamic limit 
by studying relatively small lattices. This fast convergence can be clearly seen in Fig. [4], 
where we plotted the energy per particle of a Hubbard model with U = 1 aX half filling, as a 
function of the lattice size L. This result seems to indicate that the method is more or less 
size extensive for the Hubbard model, which is surprising since in general v2DM is not size 
extensize |l5l HO] . 

Another thing to remark in Fig. |3] is that for the 20-site lattice, the difference between 
XQQ and XQQT is rather small for the half-filled lattice, but larger for the other fillings, 
and that the difference gets larger when U increases. For the 50-site lattice we see that 
the XQQ result agrees nicely with the solution of the Lieb-Wu equations. For the other 
fillings no reference data are available. There are, however, two limits of the model that 
are exactly solvable. The first limit is the rather trivial case of no interaction, i.e. U = 0, 
for which the solution has already been given in Eq. ([T]). The Hamiltonian reduces to a 
single-particle operator, which means this limit is already described correctly by including 
the X and Q conditions alone. The other exactly solvable limit is when U — ?■ +oo. In this 
limit the physics of the model decouples into two independent parts, one describing the 
spin of the system, and the other the movement of the particles (this is called spin-charge 
separation [31]). This decoupling shows up in the Bethe-ansatz wave function: the charge 
degrees of freedom are described by a Slater determinant of spinless fermions, whereas the 
spin degrees of freedom become equivalent to a spin-^ Heisenberg model. The single-particle 
energy spectrum changes slightly compared to Eq. ^ because the boundary conditions for 
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FIG. 4: Ground-state energy per particle as a function of the lattice size L, indicating how fast 
the finite size results converge to the thermodynamic limit for half filling. 

spinless fermions are periodic/antiperiodic if is even/odd PT| ITT]: 

( k = ^ if N%2 = 

ek = —2t cos k where \ , ^ ■ (40) 

|^^^(2n+lK N%2 = 1 

When the lattice is half-filled all the single-particle states are occupied, and the total energy 
sums up to zero, which is correctly described by the XQQ results in Fig. [3| Away from 



half filling, however, the energy has a finite limit which can be calculated using Eq. (40). 
From the figure we can see that the TQQ conditions do not suffice to correctly describe the 
laige-U limit. Only when the T conditions are added, the results converge to the right limit. 
Calculations at very large values of U have been performed that confirm this statement, and 
these results are shown in Table [H 



Correlation functions 



Two-particle correlation functions are important quantities in the analysis of lattice sys- 
tems, because they usually display the physics (for instance the appearance of magnetism) 
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TABLE I: Energy per site of the v2DM calculations away from half filling at large values of U, 
compared to the MPS and the Bethe-ansatz results where available. 



L 


N 


u 




TQGT 


vMPS 


exact 


L 


N 


u 


TQG 


exact 






50 


-1.2259 


-1.0804 


-1.0488 


* 






50 


-1.2272 


* 




12 


100 


-1.2177 


-1.0646 


-1.03116 


* 




30 


100 


-1.2191 


* 


20 




CXD 


* 


* 


* 


-1.0008 


50 




oo 


* 


-1.0008 




50 


-0.7972 


-0.5458 


-0.5205 


* 




50 


-0.7974 


* 




16 


100 


-0.7860 


-0.5179 


-0.49513 


* 




40 


100 


-0.7862 


* 






oo 


* 


* 


* 


-0.4639 






oo 


* 


-0.4671 



present in the system . In this Section we show that in our approach, these correlation 
functions are easily extracted from the 2DM, and compare our results to those in [301 El]- 
Charge correlation The two-particle charge correlation function is defined as: 



C{t) 



) — y~^(Q|(jQj<TQ|+r-;a'Qj+r;CT') 



(41) 



in which the notation (.) denotes the expectation value. The function is independent of the 
specific choice of the index j because of the periodic boundary conditions. The expression 



in Eq. (41) can be written in terms of the QiT) matrix: 

C(^) = ^ ^(r)j<Tj(7;(i+r)CT'(j+r)cr' , 
ucr' 

and in fact only the singlet part of the Q matrix appears: 



(42) 



C{r) = 2 g{Tf 



ji\{j+r){j+r) ■ 



(43) 



In translationally invariant systems one usually takes the Fourier transform of the correlation 
function, 

C{k) = e*'C(r) = 2J2J2 ^(H^!....., • (44) 

In Figs, [sj |6] and [7| C{k) has been plotted for ^, ^ and half filling respectively, using 
both XQQ and XQQT conditions. Comparing the XQQ with the XQQT results the same 
trends can be noticed as for the energy and the momentum distributions. For half filling 
(Fig. |7|) the XQg and XQQT results are in nice agreement. Moving away from half-filling 
(Figs. Is] and [6]) there is only agreement for small values of U . For larger values of U strange 
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FIG. 5: Two-particle charge correlation function C{k), as a function of momentum, for a ^ filled 
lattice and various values of on-site repulsion U, using IQG (top) and XQQT (bottom) conditions. 
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FIG. 6: Two-particle charge correlation function C{k), as a function of momentum, for a ^ filled 
lattice and various values of on-site repulsion U, using IQG (top) and XQQT (bottom) conditions. 



19 



oscillations appear in the XQQ results. So in this limit not only the energy, but the entire 
physical content of the IQQ-2DM cannot be trusted. This is once again an indication that 
the XQQ conditions fail to describe the strong-correlation limit away from half-filling. The 
IQGT results compare well, both in shape and magnitude, with the results from Quantum 
Monte Carlo [30], and the Bethe-ansatz results in the strong-correlation limit |31j . 
Spin correlation The two-particle spin-correlation function defined as: 

S{r) = (5^5^+'') = {a]^aj„a]_^^^..^,aj+r;a') , (45) 

era' 

can be expressed ClS db function of the Q matrix: 

S{r) = '^(r(r'g(T)j„j^.(^j+r)a'ij+r)a' ■ (46) 

Written in terms of the spin-coupled Q matrix, only the triplet S = 1 part contributes: 

Sir) = ^^(r)j,o.,.)(,^.) . (47) 
Fourier transforming S{r) yields the momentum dependent spin-correlation function: 

S{k) = Y: e^''S{r) = ^ E E ^5..... • (48) 

for ^, and half filling respectively, using 
both XQQ as XQQT conditions. Unsurprisingly, a good agreement is observed between the 
XQQ and XQQT results for small values of U . At larger values of f/, away from half filling, 
results are very poor with the XQQ conditions, especially in the :^-filled case, where the cor- 
relation function is wildly oscillating. More surprising, is that the spin-correlation function 
for the half-filled lattice in the large-f/ limit is also incorrect in the XQQ approximation. In 
the strong-correlation limit, for half-filling, the spin part of the Hubbard model is identical 
to the Heisenberg model [31], for which the spin-correlation function has a singularity at two 
times the Fermi momentum 2kF = vr, which is exactly what we see in the XQQT results. 
Below half-filling the singularity in the large-?/ limit splits and shifts to smaller values of k, 
as obserbed in the XQQT figures |8] and |9] This is in agreement with the results in [31J and 
[5U] . In conclusion we can say that the 2DM obtained with the XQQT conditions, correctly 
describes the physics that governs the spin-correlation function, whereas the XQQ conditions 
do not. It is also important to note, that even though the XQQ results for the energy are 
good for the half-filled lattice, the 2DM is flawed, because the spin-correlation function is 
not correctly described. 
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We have plotted this object in Figs, pi and 
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CONCLUSION 



In this article we have studied the one-dimensional Hubbard model at various fillings using 
the v2DM method with both two- and three-index constraints. We have shown that it is 
possible to obtain a huge reduction in the computational cost of a basic matrix computation 
by exploiting all available symmetries, i.e. spin, translation invariance and space-inversion 
parity. Using this reduction it was possible to compare the computational scaling of different 
semi definite programming algorithms with increasing lattice size. We found that, for this 
particular type of problem, the boundary point method outperforms interior point methods 
by several orders of magnitude. To gauge the quality of the variationally obtained 2DM we 
compared several ground-state properties to reference results. We found that, for half filling, 
the ground-state energy is well described by the two-index conditions. When moving away 
from half filling, however, we see that the three-index conditions are needed to obtain decent 
results. An explanation of why this happens has been the subject of a different article |35] . 
The need for three-index constraints was even more obvious when we looked at the spin and 
charge correlation functions. It was also seen that, even though the energy was relatively 
correct for the half-filled lattice, the 2DM was fiawed, because the spin correlation function 
was incorrect. This study shows that the exploitation of symmetry opens the possibility 
for a study of the two-dimensional Hubbard model for relevant lattice sizes. To obtain a 
decent accuracy, however, it will be necessary to include the three-index constraints, which 
is computationally hard. One way around this was set forward in [35J with the use of lifting 
conditions [Ml HE] . 
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FIG. 8: Two-particle spin correlation function S{k), as a function of momentum, for a filled 



lattice and various values of on-site repulsion f7, using XQQ (top) and XQQT (bottom) conditions. 
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FIG. 9: Two-particle spin correlation function S{k), as a function of momentum, for a ^ filled 
lattice and various values of on-site repulsion U, using IQG (top) and XQQT (bottom) conditions. 
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FIG. 10: Two-particle spin correlation function S{k), as a function of momentum, for a half-filled 
lattice and various values of on-site repulsion U, using ZQQ (top) and XQQT (bottom) conditions. 
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